
    ##########################################################################################
    ### R code for Figures 1 and 3 of Supplemental Material
    ### : Monte Carlo Experiment to Show that Interest Group Ratings tend to overemphasize polarization:
    ##########################################################################################
    
    ### Note: This is one-time simulation. As such, the figure might look slightly different from the figure in Supplemental Materials
    ###       Check Figure 4 for the results based on multiple simulations.
    
    library(mvtnorm)
    library(MCMCpack)

    ### I. Data Generating Process:
    #############################################################

    ### 1. Fix the numbers of legislators and roll call votes.

    M <- 100           # Number of roll call votes
    n.D <- 258
    n.R <- 176
    N <- n.D + n.R       # Number of Legislators

    G <- 9
    D.mean <- -c(0, 1, 2, 0, 1, 2, 0, 1, 2)
    R.mean <- -D.mean
    
    upper.margins <- c(0.6, 0.6, 0.6, 0.7, 0.7, 0.7, 0.8, 0.8, 0.8)
    lower.margins <- c(0.4, 0.4, 0.4, 0.3, 0.3, 0.3, 0.2, 0.2, 0.2)
 
    IRT   <- matrix(NA, N, G)
    ADA   <- matrix(NA, N, G)
    n.close.votes <- rep(NA, G)
    Vote.Margins <- matrix(NA, M, G)


    ### Number of Simulation: Vary the means the two ideal point distributions and the vote margins for ADA

    ### Generate Ideal Points and Fix Them
    Ideal.half <- matrix(NA, N, 3)
    Cutoff.half <- matrix(NA, M, 3) 
    for(g in 1:3){
    Ideal.half[,g] <- c(rnorm(n.D, mean=D.mean[g], sd=1), rnorm(n.R, mean=R.mean[g], sd=1))
    Cutoff.half[,g] <- rnorm(M, mean=0, sd=1)
    }
    Ideal <- cbind(Ideal.half, Ideal.half, Ideal.half)
    Cutoff <- cbind(Cutoff.half, Cutoff.half, Cutoff.half)
    
    for (g in 1:G){
    Ideal[c(1,N),g] <- c(-2.5, 2.5) # Ideal points for constraints
    }
      
    
    
    ### LOOP over simulations
    ##############################################################
    
    for (g in 1:G){

    ### Generate voting data
        
    Y <- matrix(NA, N, M)
    row.names(Y) <- paste("L", 1:N, sep="")
    
    # the following code is same as "Y ~ dbern(P)" in WinBUGS
    for (i in 1:N){
    for (j in 1:M){
    if (Ideal[i,g] >= Cutoff[j,g]) Y[i,j] <- 0.0
    else Y[i,j] <- 1.0
    }}

    # Compute vote margins
    for(j in 1:M){
    Vote.Margins[j,g] <- sum(Y[,j])/N
    }

    # Select votes with close margins 
    close.votes <- rep(NA, M)
    for (j in 1:M){    
    if( Vote.Margins[j,g] > lower.margins[g] & Vote.Margins[j,g] <= upper.margins[g]) close.votes[j] <- 1.0
    else close.votes[j] <- 0.0
    }

    Y.close <- Y[, close.votes==1]
    n.close.votes[g] <- n.Y.close <- ncol(Y.close)
  
    ### Estimation
    ## First, IRT estimation using MCMCirt1d
                                         
    MCMC <- MCMCirt1d(Y, theta.constraints=list(L1="-", L434="+"),  
                      burnin=500, mcmc=2000, thin=2, store.item=FALSE, store.ability=TRUE)   
    iter <- nrow(MCMC)

    # Save the mean ideal point estimates
    IRT[,g] <- apply(MCMC[,1:N], 2, FUN=mean)
  
    ## Second, ADA score 
    ADA.row  <- apply(Y.close, 1, FUN=sum)
    ADA[,g] <- 100*(ADA.row/n.Y.close)
    }
    ###############################################################
    ### End of Loop ##############################################



    #################################################
    ## Loop for producing plots


    postscript("Supplemental_Figure1.eps", height=3, width=8, horizontal=F)
    par(mfrow=c(1, 3))
    for(g in 1:3){
    plot(Ideal[,g], IRT[,g], type='p', col='black', 
         xlim=c(min(Ideal[,g]), max(Ideal[,g])), 
         ylim=c(min(IRT[,g]), max(IRT[,g])), 
         xlab="True Ideal Points", ylab="IRT Estimate")
    segments(min(Ideal[,g]), min(IRT[,g]), max(Ideal[,g]), max(IRT[,g]), lty=1)
    title(main=paste("DEM ~ N(", D.mean[g], ",1); GOP ~ N(", R.mean[g], ",1)", sep=""), font.main=1)  
    }
    dev.off()


    ADA <- -ADA
    
    postscript("Supplemental_Figure3.eps", height=8.5, width=8, horizontal=F)
    par(mfrow=c(3, 3))
    for(g in 1:G){
    plot(Ideal[,g], ADA[,g], type='p', col='black', 
         xlim=c(min(Ideal[,g]), max(Ideal[,g])), 
         ylim=c(min(ADA[,g]), max(ADA[,g])), 
         xlab="True Ideal Points", ylab="ADA Estimate")
    segments(min(Ideal[,g]), min(ADA[,g]), max(Ideal[,g]), max(ADA[,g]), lty=1)
    title(main=paste("DEM ~ N(", D.mean[g], ",1); GOP ~ N(", R.mean[g], ",1)", "\n Majority Size = ", 100*upper.margins[g], "% or less", sep=""), font.main=1)  
    }
    dev.off()
    


